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We  seek  to  predict  the  mechanical  properties  of  human  tissue  from  atomistic  simulations  and 
microrheological  experiments.  The  most  prevalent  protein  in  the  body  is  also  one  that  determines  much  of 
its  mechanical  properties:  collagen.  Collagen  is  a  multiscale  material,  so  requires  multiscale  modeling  for 
predictions.  It  consists  of  long  proteins  that  form  a  triple  helix.  The  triple  helices  stack  to  form  fibrils  of  a 
few  hundred  nanometers  in  diameter.  In  turn,  these  fibrils  then  bind  to  form  fibers  that  have  the  important 
mechanical  strength.  We  show  that  the  modulus  of  the  triple  helix  and  fibrils  are  very  similar,  but  that  the 
fibers  have  much  lower  modulus.  Hence,  we  postulate  that  the  inter-fibril  cross-linking  is  significantly  softer 
than  collagen.  Details  are  in  the  submitted  manuscript  Nanomechanics  of  type  I  collagen  fibrils,  and  only  a 
brief  summary  is  given  below. 

Our  overall  strategy  is  the  following  steps.  Co-principal  investigator  Joseph  Orgel  has  performed  x-ray 
diffraction  to  determine  the  atomistic  structure  of  the  stacking  in  the  fibrils.  These  structures  are  then  fed 
into  an  atomistic  molecular  dynamics  simulation.  Many  of  the  force  fields  between  atoms  are  determined 
from  ab  initio  quantum  mechanical  calculations,  and  these  are  then  used  to  determine  the  stresses  in  the 
fibril  material  as  function  of  strain.  The  simulations  use  periodic  boundary  conditions  to  assume  that  the 
fibrils  are  infinitely  large.  The  results  of  the  simulation  are  used  to  determine  parameters  in  a  linear  elastic 
rod  model.  These  predictions  are  compared  to  single-fiber  mechanical  measurements  performed  using 
optical  tweezers  (see  manuscript  "Measurement  of  elastic  modulus  of  collagen  Type  I  single  fiber").  Finally, 
microrheology  is  performed  on  gels  made  of  collagen.  Most  details  are  given  in  previous  reports,  published 
papers  and  proceedings,  all  of  which  are  included  elsewhere  in  the  report.  In  the  following  sections  we 
focus  here  on  what  is  new  and  not  in  those  other  publications. 

We  also  examine  the  properties  of  oriented  and  unoriented  gels.  We  have  developed  techniques  to 
quantify  the  anisotropy  of  the  gel,  whose  structure  was  created  by  periodic  deformations  during 
crosslinking.  We  also  show  that  two-point  microrheology  can  be  used  to  measure  the  modulus  of  gels 
whose  microstructure  is  larger  than  bead  size.  These  details  are  given  below. 


Molecular  dynamics  simulation 

Collagen  is  the  primary  constituent  of  animal  connective  tissue.  A  defining  feature  of  fibrillar  collagen  is  its 
axial  periodicity  visible  in  transmission  electron  microscopy  as  alternating  dark  and  light  bands.  The 
repeating  unit,  D,  is  67  nm  long,  with  the  dark  band  making  up  54%  of  the  repeating  length.  This  periodicity 
of  dark/light  bands  reflects  an  underlying  packing  of  constituent  triple-helix  polypeptide  monomers  wherein 
the  dark  bands  represent  gaps  between  axially  adjacent  monomers.  This  organization  of  gap/overlap 
regions  is  present  in  the  microfibrillar  model  of  collagen  obtained  from  fiber  diffraction  at  55  Angstrom 
resolution. 

Recent  molecular  dynamics  simulations  of  this  model  under  zero-stress  conditions,  however,  predict  that 
the  D-band  shrinks  by  19%  [1].  Consequently,  the  mechanical  properties  obtained  from  this  study 
correspond  to  those  of  distorted  microfibrils.  Here  we  evaluate  systematically  the  effect  of  several  physical 
parameters  on  D-band  shrinking.  Using  the  force  field  employed  in  the  reported  study,  we  find  that, 
irrespective  of  the  temperature/pressure  coupling  algorithms,  assumed  salt  concentration  or  hydration 
level,  and  whether  or  not  the  monomers  are  cross-linked,  the  D-band  shrinks  considerably.  This  shrinkage 
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results  not  from  the  sliding  of  axially  parallel  monomers  across  one  another,  but  rather  it  is  associated  with 
bending  and  widening  of  individual  monomers.  The  distributions  of  backbone  dihedrals  are  also  different 
from  those  in  the  fiber  diffraction  model. 

Employing  a  force  field  whose  backbone  dihedral  energy  landscape  matches  more  closely  with  CCSD(T) 
theory  results  in  a  small  D-band  shrinking  of  <3%.  Since  this  force  field  is  also  known  to  perform  better 
against  other  experimental  structural  data,  it  appears  that  the  large  shrinkage  observed  in  earlier 
simulations  is  a  force-field  artifact.  However,  certain  atomic-level  details,  like  glycosylation  sites,  which  are 
indiscernible  from  the  electron  densities  of  collagen  and  collagen-related  peptides,  still  remain  to  be 
incorporated  into  the  atomically  detailed  model  of  collagen.  It  is,  therefore,  plausible  that  their  inclusion 
might  prevent  large  force-field-dependent  variations  in  D-band  length.  Consequently,  we  cannot  yet 
completely  rule  out  D-band  shrinking  to  be  an  artifact  of  model  construction. 

Fibrils  of  type  I  collagen  assemble  from  polypeptide  triple  helices.  Higher-order  structures,  such  as  fibril- 
bundles  and  fibers,  assemble  from  fibrils  in  the  presence  of  other  collagenous  and  non-collagenous 
molecules.  A  compilation  of  results  from  several  different  experiments  indicates  that  the  Young's  moduli  of 
fibrils/fibril-bundles  are  an  order  in  magnitude  smaller  than  those  of  triple  helices,  indicating  that 
fibrils/fiber-bundles  are  less  resistant  to  axial  deformation  compared  to  triple  helices.  To  understand  the 
molecular  basis  for  this  difference,  we  carry  out  all-atom  molecular  dynamics  simulations  of  the  smallest 
repeating  crystallographic  unit  of  a  fibril  under  periodic  boundary  conditions.  This  system  is  structurally 
equivalent  to  an  infinite  microfibril.  Assuming  that  linear  response  theory  applies  to  axial  strains  smaller 
than  5%,  we  find  that  the  stress-strain  relationship  yields  a  Young's  modulus  of  2.34  GPa,  which  is  within  the 
range  estimated  for  triple  helices.  This  similar  magnitude  in  modulus  suggests  that  the  smaller  resistance  of 
collagen  fibrils/fiber-bundles  to  axial  deformation  does  not  stem  from  the  packing  of  polypeptide  triple 
helices  in  a  microfibril.  Instead,  it  stems  from  the  assembly  of  microfibrils  into  finite-sized  fibrils  and/or  the 
assembly  of  fibrils  into  fibril-bundles. 


Anisotropic  mechanical  response  of  oriented  collagen  gels 

Tissue  and  tissue  scaffolding  made  from  collagen  is  usually  oriented,  rather  than  isotropic.  The  cells 
appear  to  be  sensitive  to  their  anisotropic  mechanical  properties,  which  we  seek  here  to  quantify.  Towards 
this  end,  we  use  microrheology  to  examine  the  anisotropic  fluctuations  of  Brownian  beads  embedded  in 
networks  that  are  formed  under  flow.  To  impose  the  flow,  we  designed  and  built  a  small  device  that  could 
induce  flow  inside  our  commercial  device  for  bead  tracking  and  manipulation. 

The  JPK  NanoTracker  [2]  was  used  to  record  thermal  fluctuations  of  2  pm  polystyrene  beads 
(Polybead  Polystyrene  Microspheres,  mean  diameter  2.0  pm,  standard  deviation  ±5  %,  Polysciences,  Inc., 
Warrington,  PA),  embedded  into  1  mg/ml  self-assembled  rat  tale  type  I  collagen  network  (Corning® 
Collagen  I,  Rat  Tail,  Cornong  Inc.,  Corning,  NY) 

To  create  anisotropy  in  the  sample  we  used  a  custom-built  device  to  produce  a  large  amplitude 
oscillatory  shear  (LAOS)  flow  in  the  sample  chamber  during  the  gelation  process  (Figure  1).  The  sample  (7-8 
pi)  was  placed  between  microscope  coverslips  (18x18  mm  and  24x55  mm,  Nol,  Thermo  Fisher  Scientific, 
Fermont,  CA)  (Figure  2).  The  top  coverslip  was  repeatedly  moved  by  a  linear  motor  (PiezoMotor  Piezo  LEGS® 
Linear  6N)  with  an  amplitude  of  100  pm  and  maximum  speed  10  pm/s  during  the  first  3  hours  of 
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fibrogenesis.  To  avoid  a  complicated  secondary  flow  pattern  during  the  shear  deformation,  only  a  small 


rectangular  portion  of  the  sample  chamber  was  actually  filled  with  the  sample  (Figure  3). 

Figure  1  Render  of  shear  chamber. 

Figure  2  Shear  chamber  cross-section 
Figure  3  Shear  chamber  top  view  (Sizes  in  mm) 

Microrheological  experiments  were  performed  as  follows.  The  collagen  gel  was  injected  into  the 
sample  chamber  and  left  for  3  hours  under  LAOS  flow.  An  optical  trap  was  focused  on  the  desired  bead  and 
calibrated.  After  that,  data  of  the  bead  position  was  collected  simultaneously  for  X  and  Y  channels  at  the 
rate  of  100  kFIz  for  60  seconds.  Calibration  and  data  collection  steps  were  repeated  for  5-10  different  beads. 
The  sample  chamber  was  turned  90°  counterclockwise  and  the  whole  process  of  the  calibration  and  the 
data  collection  was  repeated  for  exactly  the  same  beads.  The  last  step  is  important  for  a  reproducibility 
check  and  for  evenly  distributing  all  the  internal  and  external  noise  sources  among  X  and  Y  channels. 


Results  and  Discussion 
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There  are  two  competitive  processes  occurring  in  the  sample  during  LAOS.  One  process  is  the 
orientation  of  fibers  along  the  direction  of  the  shear  as  described  by  [3],  The  other  process  caused  by  LAOS 
is  weakening  of  the  network  in  the  flow  direction  [4],  As  described  in  [4],  the  weakening  occurs  in  a  mature 
network  due  to  destruction  of  cross-linking  bounds  between  fibers.  The  effect  was  shown  for  relatively 
small  oscillatory  strains  (0.02  -  0.2)  and  large  step  strain  (0.7)  followed  by  long  relaxation  time  (tens  of 
minutes).  However,  during  fibrogenesis  we  noticed  that  cross-links  are  easily  reformed  due  to  the  presence 
of  proteoglycans  and  therefore  the  fibers  are  just  orient  in  the  direction  of  LAOS  flow  (direction  along  y 


axis  on  Figure  4,  b). 

Figure  4  Unoriented  (a)  and  oriented  (b)  collagen  networks.  For  the  oriented  network  (b),  the  x  axis 
is  vorticity  direction  and  the  y  axis  is  flow  direction. 

Observing  fibrogenesis  under  LAOS  in  real  time,  we  were  able  to  see  how  newly  formed  fibers 
become  oriented  in  the  flow  direction  ( y  on  Figure  4,  b)  and  then  bonded  to  the  existing  network. 
Moreover,  fibers  aggregate  in  some  regions  leaving  other  regions  empty.  For  example,  the  left  side  of  Figure 
4,  b  is  empty  of  fibers  while  the  middle  and  right  sides  contain  oriented  fiber  bundles  which  are  visibly  more 
dense  than  the  unoriented  network  Figure  4,  a. 

The  mean-square  displacement  (MSD)  of  the  particle  positions  on  two  perpendicular  axes  (x  and  y. 
Figure  4)  was  calculated  for  8  different  particles  embedded  in  unoriented  (Figure  5  and  Figure  6)  and 
oriented  (Figure  7  and  Figure  8)  collagen  networks. 

One  notices  significant  variance  between  the  MSDs  for  different  particles  for  both  x  and  y  channels. 
Such  variance  was  reported  before  [5]  and  is  probably  caused  by  the  slight  inhomogeneity  within  the 
formed  fiber  bundle  (right  side  of  Figure  4,  b).  To  avoid  inhomogeneity  one  can  try  increasing  the  collagen 
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concentration.  However,  higher  collagen  concentration  will  lead  to  significant  stiffening  of  the  network, 
which  in  turn  decreases  the  amplitude  of  the  particle  motion,  possibly  to  the  resolution  limit.  The  latter 
might  be  especially  important  for  the  oriented  networks,  which  tend  to  already  have  higher  local  collagen 
concentrations  due  to  the  formation  of  aggregates  (bundles). 

Comparing  MSD  plots  for  x  and  y  channels,  one  sees  that,  despite  the  variation  between  different 
particles,  the  plots  look  similar.  To  further  analyze  the  difference  between  x  and  y  displacements,  we  plot 
their  relation  in  Figure  9.  For  an  ideal  isotropic  network  MSD  x/MSD  y  should  always  equal  1.  However,  even 
for  our  unoriented  1  mg/ml  collagen  gel  there  are  some  deviations  from  1.  That  could  be  explained  by  the 

MSDx 

small  concentration  of  fibers  in  the  gel.  Another  feature  of  Figure  9  is  that  the  quantity - , 

MSD  y 

which  can  be  explained  increased  contribution  of  inertia  of  the  particle  and  the  buffer  at  high  frequencies. 

To  quantify  the  anisotropy  of  the  sample,  we  calculate  the  average  and  the  standard  deviation  of 
MSD  x/MSD  y  for  different  particles  at  a  lag  time  of  0.1  s  which  turned  out  to  be  1.49±0.83. 

Analyzing  MSD  plots  for  the  oriented  network  (Figure  7  and  Figure  8),  one  sees  that  MSDs  for  the 
oriented  network  is  generally  smaller  than  for  unoriented.  This  is  consistent  with  the  oriented  network 
having  higher  local  collagen  concentrations,  as  was  discussed  earlier. 

Another  feature  is  that  plots  of  MSD  y  (which  correspond  to  fluctuation  of  particles  along  the 
preferred  fiber  orientation  in  bundle  on  Figure  4,  b)  for  different  particles  are  generally  below  plots  of  MSD 
X  for  same  particles  (which  correspond  to  fluctuation  of  particles  perpendicular  to  the  preferred  fiber 
orientation  in  the  bundle  on  Figure  4,  b).  This  appears  to  be  an  effect  of  anisotropy  in  the  fiber  bundles  as 
motion  is  more  restricted  in  the  direction  y  on  Figure  4. 

Lastly,  a  significant  amount  of  noise  appears  in  both  Figure  7  and  Figure  8.  This  noise  has  a  distinct 
spectral  signature  and  is  believed  to  be  a  result  of  external  vibrations.  We  took  serious  efforts  in  reducing 
the  noise  propagated  to  the  samples,  but  we  were  only  able  to  reduce  it  below  10'^^  m^. 

To  quantify  the  anisotropy  of  the  oriented  sample,  we  again  calculate  the  average  and  the  standard 
deviation  of  MSD  x/MSD  y  for  different  particles  at  the  lag  time  of  0.1  s  which  turned  out  to  be  3.33±1.26,  a 
clear  quantification  of  anisotropy  in  the  sample. 


5 


Figure  5  Mean-square  displacement  along  the  x  axis  of  the  chamber  for  8  different  particles  (P1-P8) 


embedded  into  non-oriented  collagen  network. 


Figure  6  Mean-square  displacement  along  the  y  axis  of  the  chamber  for  8  different  particles  (P1-P8) 
embedded  into  non-oriented  collagen  network. 
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Figure  7  Mean-square  displacement  along  the  x  axis  of  the  chamber  for  8  different  particles  (P1-P8) 
embedded  into  oriented  collagen  network. 


Figure  8  Mean-square  displacement  along  the  y  axis  of  the  chamber  for  8  different  particles  (P1-P8) 
embedded  into  oriented  collagen  network. 
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Figure  9  Relation  between  MSD  for  x  and  y  displacements  (MSD  x  /MSD  y)  for  8  different  particles 


(P1-P8)  embedded  into  non-oriented  collagen  network. 

Figure  10  Relation  between  MSD  for  x  and  y  displacements  (MSD  x  /MSD  y)  for  8  different  particles 
(P1-P8)  embedded  into  oriented  collagen  network. 


Two-bead  microrheology 

We  have  also  developed  detailed  methods  for  analyzing  the  results  of  two-bead  microrheology[6].  Roughly 
speaking,  one  measures  the  cross-correlations  of  the  Brownian  motion  of  two  beads  that  interact 
hydrodynamically  in  a  viscoelastic  medium.  In  single-point  microrheology,  a  single  bead  is  followed,  but  the 
bead  must  be  larger  than  the  microstructure  of  the  medium.  If  the  bead  is  too  large,  it  is  not  sufficiently 
subject  to  Brownian  motion,  and  no  signal  can  be  detected.  Two-point  microrheology  overcomes  these 
limitations  by  using  beads  that  are  small  enough  to  experience  Brownian  forces,  but  significantly  separated 
to  sample  the  microstructure  as  a  continuum. 
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Prior  analysis  neglected  two  effects:  the  finite  time  for  waves  to  travel  between  the  two  particles,  and  the 
fact  that  a  single  wave  can  bounce  between  the  two  particles  more  than  once.  We  have  included  these 
effects  and  designed  an  analysis  that  is  no  more  complicated  than  the  traditional  approach  [7],  [8],  We 


apply  those  analyses  here  to  our  own  data  in  a  prototype  polymer  solution,  as  shown  in  Figure  11. 

Figure  11  Mean-squared  displacement,  and  cross-correlations  in  two-bead  microrheology.  The 
material  is  2%  PAM(polyacrylamide)  dissolved  in  water.  The  beads  have  radius  i?  =  0.765//m  and  are 
separated  by  AR  .  The  solid  lines  are  a  fit  used  to  convert  these  data  from  the  time  domain  to  the  frequency 
domain. 

These  data  are  then  used  to  estimate  the  dynamic  modulus  of  the  polymer  solution  and  compared  to  the 
same  quantity  measured  in  a  mechanical  rheometer.  The  results  are  shown  in  Figure  12. 
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Figure  12  Dynamic  modulus  measured  by  a  mechanical  rheometer  in  parallel  plates  (Macrorheology),  and 
that  measured  using  single  and  two-point  microbead  rheology. 


Semiflexible  fiber  mechanics 

Further  coarse  graining  requires  statistical  mechanics.  Flere  we  model  the  fibers  as  fluctuating  elastic  fibers 
with  uniform  Young's  modulus,  the  so-called  "wormlike  chain"  [9].  Great  strides  were  recently  made  in 
these  derivations  and  calculations  by  Spakowitz  and  coworkers  [10].  In  order  to  build  mathematical  model 
of  networks  of  these  strands,  we  require  the  free  energy  of  a  strand  which  interacts  with  other  strands 
through  topological  constraints  and  cross-links.  Between  these  constraints,  the  strands  are  relatively  free  to 
fluctuate  from  Brownian  forces.  Flowever,  the  strand  orientations  at  these  constraints  must  be  continuous, 
which  couples  the  free  energy  of  two  adjacent  chain  segments. 

Using  the  approach  taken  by  Koslover  and  Spakowitz  [11],  we  found  the  Green's  function  (equivalent  to  the 
free  energy)  of  two  strands  with  fixed  length  ( L  ),  fixed  end-to-end  vector  ( Q  ),  and  prescribed  orientation 
at  their  junction  ( u ) 


An 


(1) 


where  the  Y  are  spherical  harmonics,  the  j  are  spherical  Bessel  functions,  and  the  G°’°  are  continued 
fractions. 

As  an  example,  we  consider  the  buckling  of  this  strand  compared  to  the  prediction  of  Euler  buckling,  which 
does  not  include  fluctuations.  Schematic  representation  of  the  buckled  chain  and  notations  are  shown  in 
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Figure  13.  Figure  14  shows  the  critical  force  required  for  buckling  relative  to  the  Euler  prediction.  The  result 
shows  that  fluctuations  actually  stabilize  the  fiber.  It  appears  that  for  short  chains,  the  critical  force  scales 
linearly. 


U2lp 

Figure  13  Schematic  representation  of  the  chain  modeled  as  2  equal  strands. 

Figure  14  Critical  buckling  force  vs  contour  length  of  the  chain 

We  also  looked  at  the  distribution  of  the  tangent  orientation  u  at  the  middle  of  the  buckled  strand. 
Examples  of  this  distribution  are  shown  in  Figure  15,  Figure  16  and  Figure  17. 

For  10%  compression  (Rl L  =  0.9  )  the  chain  is  almost  straight  and  the  distribution  is  sharply  peaked  at  the 
direction  of  the  end-to-end  vector,  allowing  only  small  isotropic  fluctuation  around  it  (Figure  15).  For  50% 
compression  the  orientation  of  the  middle  point  goes  out  of  the  plain,  indicating  that  chain  configurations 
can  become  helical  (Figure  16).  This  suggests  that  2D  models  would  not  be  sufficient  to  describe  buckling  of 
semiflexible  fibers,  and  shows  the  importance  of  fluctuations. 
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Figure  15  Distribution  of  orientation  angle  for  the  middle  point  on  the  chain  L  !  21^=2  , 


R/ L  =  0.9  ( I L  =  0.06  ).  Color  map  is  the  same  for  all  a 

Figure  16  Distribution  of  orientation  angle  for  the  middle  point  on  the  chain  L  !  21^  =  2 , 
R! L  =  Q. 5  ( I L  =  0.3  ).  Color  map  is  the  same  for  all  a  . 


Figure  17  Distribution  of  orientation  angle  for  the  middle  point  on  the  chain  L  !  21^  =  2 , 
R! L  =  Q.\  ( I L  =  0.33  ).  Color  map  is  the  same  for  all  a  . 


The  free  energies  are  then  suitable  for  building  multi-chain  or  single-chain  mean-field  models  of  entangled 
cross-linked  semi-flexible  networks,  such  as  those  that  make  up  tissue. 
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